* Nybom and Stuhler, Interpreting Trends in Intergenerational Mobility - Figure A1
* Simulate mobility trend for cohorts vs. generations example

* prepare Stata
clear all
set more off 
set obs 1000
cd "/Users/janstuhler/Dropbox/MartinJanProjects/Trend_Project/Submissions/JPE/Final version/Simulations/Figure A1"

* import distribution of father's age at birth from Swedish data at year 1960 (drop fathers with f_ageatbirth<18 and ... >50, about 1% of cases)
scalar dens_18 = 0.0037 
scalar dens_19 = 0.0078 
scalar dens_20 = 0.0145 
scalar dens_21 = 0.0230 
scalar dens_22 = 0.0333 
scalar dens_23 = 0.0390 
scalar dens_24 = 0.0452 
scalar dens_25 = 0.0492 
scalar dens_26 = 0.0550 
scalar dens_27 = 0.0560 
scalar dens_28 = 0.0579 
scalar dens_29 = 0.0593 
scalar dens_30 = 0.0568 
scalar dens_31 = 0.0528 
scalar dens_32 = 0.0511 
scalar dens_33 = 0.0481 
scalar dens_34 = 0.0444 
scalar dens_35 = 0.0426 
scalar dens_36 = 0.0389 
scalar dens_37 = 0.0341 
scalar dens_38 = 0.0311 
scalar dens_39 = 0.0293 
scalar dens_40 = 0.0268 
scalar dens_41 = 0.0192 
scalar dens_42 = 0.0166 
scalar dens_43 = 0.0144 
scalar dens_44 = 0.0117 
scalar dens_45 = 0.0096 
scalar dens_46 = 0.0084 
scalar dens_47 = 0.0067 
scalar dens_48 = 0.0059 
scalar dens_49 = 0.0041 
scalar dens_50 = 0.0034 

* Density young and old parents
scalar dens_young=0
forval i=18/30 {
  scalar dens_young=dens_young+dens_`i'
}
scalar dens_old=0
forval i=31/50 {
  scalar dens_old=dens_old+dens_`i'
}
di dens_young
di dens_old
	
	
* set parameter values before and after structural change
scalar rho1 = 0.2
scalar lambda1 = 0.8
scalar gamma1 = 0.3
scalar varu=0.844736842	

* pre-shock
forvalues c=1700(1)1959 {
	scalar rho_`c' = 0.2
	scalar lambda_`c' = 0.8
	scalar gamma_`c' = 0.3
  
}
* post-shock
forvalues c=1960(1)2080 {
	scalar rho_`c' = 0.5
	scalar lambda_`c' = 0.8
	scalar gamma_`c' = 0.3	
}
* smooth shock?
/*scalar rho_1960 = 0.23
scalar rho_1961 = 0.26
scalar rho_1962 = 0.29
scalar rho_1963 = 0.32
scalar rho_1964 = 0.35
scalar rho_1965 = 0.38
scalar rho_1966 = 0.41
scalar rho_1967 = 0.44
scalar rho_1968 = 0.47
scalar rho_1969 = 0.50*/

** Method (1): Start from steady state, then recursive solution
gen year=.
gen beta=.
gen cov_ey=.
gen cov_ey_p=.
gen var_y=.
gen var_y_p=.

gen beta_young=.
gen cov_ey_p_young=.
gen var_y_p_young=.

gen beta_old=.
gen cov_ey_p_old=.
gen var_y_p_old=.

* Steady state covariances
local obs = 0
forvalues c=1860(1)1959 {

	local obs = `obs'+1
	replace year=`c' in `obs'
	replace cov_ey = rho1 / (1-gamma1*lambda1) in `obs' 
	replace var_y = 1 in `obs' 	
}
* Recursive covariance and variances
forvalues c=1960(1)2080 {	

	local obs = `obs'+1
	replace year=`c' in `obs'
	
	replace cov_ey = 0 in `obs'
	replace var_y = 0 in `obs'
	
	forvalues age=18(1)50 {
		local c1 = `c'-`age'
		replace cov_ey = cov_ey + dens_`age'*(gamma_`c'*lambda_`c'*cov_ey[_n-`age'] + rho_`c')   in `obs'
		replace var_y = var_y + dens_`age'*(gamma_`c'^2*var_y[_n-`age']+ rho_`c'^2*1 + 2*gamma_`c'*rho_`c'*lambda_`c'*cov_ey[_n-`age']+varu)  in `obs'	// Law of total variance, assuming/exploiting that E[Y|parental age]=0)
	}	
}
* Elasticities
local obs = 0
forvalues c=1860(1)2080 {
	
	local obs = `obs'+1
	if `c'>=1910 {
		* IGE
		replace cov_ey_p = 0  in `obs'
		replace var_y_p = 0 in `obs'
		forvalues age=18(1)50 {
			replace cov_ey_p = cov_ey_p + rho_`c'*lambda_`c'*dens_`age'*cov_ey[_n-`age'] in `obs'
			replace var_y_p = var_y_p + dens_`age'*var_y[_n-`age'] in `obs'
		}
		replace beta = cov_ey_p/var_y_p + gamma_`c'	in `obs'
		* IGE young parents
		replace cov_ey_p_young = 0  in `obs'
		replace var_y_p_young = 0 in `obs'
		forvalues age=18(1)30 {
			replace cov_ey_p_young = cov_ey_p_young + rho_`c'*lambda_`c'*(dens_`age'/dens_young)*cov_ey[_n-`age'] in `obs'
			replace var_y_p_young = var_y_p_young + (dens_`age'/dens_young)*var_y[_n-`age'] in `obs'			
		}
		replace beta_young = cov_ey_p_young/var_y_p_young + gamma_`c'	in `obs'
		* IGE old parents
		replace cov_ey_p_old = 0  in `obs'
		replace var_y_p_old = 0 in `obs'
		forvalues age=31(1)50 {
			replace cov_ey_p_old = cov_ey_p_old + rho_`c'*lambda_`c'*(dens_`age'/dens_old)*cov_ey[_n-`age'] in `obs'
			replace var_y_p_old = var_y_p_old + (dens_`age'/dens_old)*var_y[_n-`age'] in `obs'			
		}
		replace beta_old = cov_ey_p_old/var_y_p_old + gamma_`c'	in `obs'
	}
}

 
/*
** Method (2): full recursive solution (takes long to run)
gen beta_rec=.

local obs = 0
forvalues c=1900(1)2080 {

	local obs = `obs'+1

	forvalues g=1(1)5 {	
		scalar gen`g' = 0
	}
	
	* generation 1 
	forvalues age1=18(1)50 {
		local c1 = `c'-`age1'
		scalar gen1 = gen1 + rho_`c'*lambda_`c'*dens_`age1'*rho_`c1'
		
		* generation 2 
		forvalues age2=18(1)50 {
			local c2 = `c'-`age1'-`age2'
			scalar gen2 = gen2 + rho_`c'*lambda_`c'*dens_`age1'* ///
				dens_`c1'_`age2'*gamma_`c1'*lambda_`c1'* ///
				rho_`c2'

			* generation 3 
			forvalues age3=18(1)50 {
				local c3 = `c'-`age1'-`age2'-`age3'
				scalar gen3 = gen3 + rho_`c'*lambda_`c'*dens_`age1'* /// 
					dens_`c1'_`age2'*gamma_`c1'*lambda_`c1'* ///
					dens_`c2'_`age3'*gamma_`c2'*lambda_`c2'* ///
					rho_`c3'
					
				* generation 4 
				*forvalues age4=18(1)50 {
				*	local c4 = `c'-`age1'-`age2'-`age3'-`age4'
				*	scalar gen4 = gen4 + rho_`c'*lambda_`c'*dens_`age1'* /// 
				*		dens_`c1'_`age2'*gamma_`c1'*lambda_`c1'* ///
				*		dens_`c2'_`age3'*gamma_`c2'*lambda_`c2'* ///
				*		dens_`c3'_`age4'*gamma_`c3'*lambda_`c3'* ///
				*		rho_`c4'
				*}
			}
		}
	}
		
	replace beta_rec= gamma_`c' + gen1 + gen2 + gen3 + gen4 in `obs' 
}
*/

* graph 
twoway line beta_young beta_old beta year if year>=1950 & year<=2040 , ///
	xlabel(1950(10)2040) ylabel(0.3(0.05)0.52) lpattern(dash shortdash solid) ///
	scheme(stcolor_alt) legend(cols(3) order(3 1 2) lab(1 "IGE, younger parents") lab(2 "IGE, older parents") lab(3 "IGE") )
graph export "Simluation_CohortsvsGenerations_Rho_SmoothShock_2023.pdf" , replace
